Making static maps
Matt Williamson: JF edit
10/20/2021
Your assignment
The goals for this assignment are to:
- practice making maps with the various packages
- think about composition and the grammar of graphics
- critique various maps
By the end of the assignment you should have several static maps displaying the datasets we’ve used in the last few weeks.
Task 1: Show me your ugly maps!!
Find 2 examples of maps that you think are ‘bad’.
Question 1 Why are they bad? What might improve them? The first map is trying to do way too much. My interpretation is that it has taken three datasets, high school graduates, college graduates, and median household income. The infographic then combines the colors for each county. While looking at the relationship of educational attainment and income has merit as a question, combining the data in this way was not an effective way to show that relationship. It is also not clear that the county level was any more interesting than say the state level. One option to improve could have been creating a ratio of education to income and normalizing that with one color scheme so places could more readily be compared to one another.
Question 2 Rely on the Healy and Wilke texts to provide some structure to your answers. The second map (as well as the first) goes against the rule of proportional ink, the symbols are giant and its just one big blob of grey. You don’t get any information from this other than a takeaway that Belgium has a lot of churches, which isn’t particularly interesting.
Task 2: Load your libraries and build a dataframe
You can choose whichever datasets you’d like from the past several months as the subject for your mapping. You’ll need to use at least one tabular join, one spatial join, and one extraction to create the dataframe. Load the packages, the data, and make sure everything is projected here. Give me a sense for what you are hoping to map.
## Reading layer `regionalPAs1' from data source
## `/opt/data/session04/regionalPAs1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 224 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -1714157 ymin: 1029431 xmax: -621234.8 ymax: 3043412
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic_USGS_version## Reading layer `reg_pas' from data source `/opt/data/session16/reg_pas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 544 features and 32 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2034299 ymin: 1004785 xmax: -530474 ymax: 3059862
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic_USGS_version## Warning in colnames(pas.proc)[c(1, 6, 8, 10, 12, 22, 25)] <-
## colnames(pas.desig): number of items to replace is not a multiple of replacement
## lengthBring in a map of Montana and match our projections
mt <- tigris::states(cb=TRUE) %>%
filter(STUSPS == "MT") ##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|=============== | 21%
|
|================ | 23%
|
|================= | 25%
|
|=================== | 27%
|
|====================== | 31%
|
|======================= | 33%
|
|========================= | 35%
|
|========================= | 36%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================= | 41%
|
|============================== | 43%
|
|================================= | 47%
|
|=================================== | 50%
|
|=================================== | 51%
|
|===================================== | 53%
|
|======================================================================| 100%# make spatvectors
st_crs(mammal.rich)$proj4string## [1] "+proj=aea +lat_0=37.5 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"pa.vect <- as(pas, "SpatVector")
mt.vect <- as(mt, "SpatVector")
# project to match mammal richness
pa.vect <- project(pa.vect, mammal.rich)
mt.vect <- project(mt.vect, mammal.rich)
land.val.proj <- project(landval, mammal.rich)
hmi.proj <- project(hmi, mammal.rich)
#crop to match montana
mam.rich.crop <- crop(mammal.rich, mt.vect)
mt.val.crop <- crop(land.val.proj, mt.vect)
hmi.crop <- crop(hmi.proj, mt.vect)
plot(mam.rich.crop)
plot(mt.vect, add=TRUE)
plot(pa.vect, add=TRUE)
Add some tidy census variables
mt.income <- tidycensus:: get_acs(geography = "county",
variables = c(medianincome = "B19013_001", # I tried pretty hard to get this to be % poverty but kept getting errors :(
pop = "B01003_001"),
state = c("MT"),
year = 2018,
key = key,
geometry = TRUE) %>%
st_transform(., crs(mammal.rich)) %>%
select(-moe) %>%
spread(variable, estimate)## Getting data from the 2014-2018 5-year ACS## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.##
|
| | 0%
|
|========== | 14%
|
|=================== | 28%
|
|============================= | 42%
|
|======================================= | 55%
|
|================================================= | 69%
|
|========================================================== | 83%
|
|==================================================================== | 97%
|
|======================================================================| 100%let’s see if I can get this to work with poverty data…
mt.poverty <- tidycensus:: get_acs(geography = "county",
variables = c(family_poverty = "B17010_002",
pop = "B01003_001"),
state = c("MT", "WY"),
year = 2018,
key = key,
geometry = TRUE) %>%
st_transform(., crs(mammal.rich)) %>%
select(-moe) %>%
spread(variable, estimate)## Getting data from the 2014-2018 5-year ACS## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.reg.poverty <- st_make_valid(mt.poverty)
pa.summary.pov<- st_join(st_as_sf(pa.vect), mt.poverty, join = st_overlaps)
pa.summary.pov <- pa.summary.pov %>%
group_by(Unit_Nm) %>%
summarize(., family_poverty = mean(family_poverty, na.rm=TRUE),
meanpop = mean(pop, na.rm=TRUE))
#double check to see that I got the right number of rows
nrow(pa.summary.pov) ==length(unique(pas$Unit_Nm))## [1] TRUE## [1] TRUERun zonal statistics
pa.zones <- terra::rasterize(pa.vect, mam.rich.crop, field = "Unit_Nm")
mammal.zones <- terra::zonal(mam.rich.crop, pa.zones, fun = "mean", na.rm=TRUE)
hmi.zones <- terra::zonal(hmi.crop, pa.zones, fun = "mean", na.rm=TRUE)
landval.zones <- terra::zonal(mt.val.crop, pa.zones, fun = "mean", na.rm=TRUE)
#Note that there is one few zone than we have in our PA dataset. This is because we have an overlapping jurisdicition; we'll ingnore that now but it's a common problemen with using the PADUS
summary.df.pov <- pa.summary.pov %>%
left_join(., mammal.zones) %>%
left_join(., landval.zones) %>%
left_join(., hmi.zones)## Joining, by = "Unit_Nm"
## Joining, by = "Unit_Nm"
## Joining, by = "Unit_Nm"Extract and put everything together into one df
pa.summary.income <- st_join(st_as_sf(pa.vect), mt.income, join = st_overlaps)
pa.summary.income <- pa.summary.income %>%
group_by(Unit_Nm) %>%
summarize(., meaninc = mean(medianincome, na.rm=TRUE),
meanpop = mean(pop, na.rm=TRUE))
#double check to see that I got the right number of rows
nrow(pa.summary.income) ==length(unique(pas$Unit_Nm))## [1] TRUE## [1] TRUE
pa.zones <- terra::rasterize(pa.vect, mam.rich.crop, field = "Unit_Nm")
mammal.zones <- terra::zonal(mam.rich.crop, pa.zones, fun = "mean", na.rm=TRUE)
hmi.zones <- terra::zonal(hmi.crop, pa.zones, fun = "mean", na.rm=TRUE)
landval.zones <- terra::zonal(mt.val.crop, pa.zones, fun = "mean", na.rm=TRUE)
#Note that there is one few zone than we have in our PA dataset. This is because we have an overlapping jurisdicition; we'll ingnore that now but it's a common problement with using the PADUS
summary.df.inc <- pa.summary.income %>%
left_join(., mammal.zones) %>%
left_join(., landval.zones) %>%
left_join(., hmi.zones)## Joining, by = "Unit_Nm"
## Joining, by = "Unit_Nm"
## Joining, by = "Unit_Nm"colnames(summary.df.inc)## [1] "Unit_Nm" "meaninc" "meanpop" "geometry"
## [5] "Value" "landvalrastmtwy" "hm_fsum3_270"Task 3: Build a map with Tmap
Practice making a quick map with tmap.
tm_shape(mam.rich.crop) +
tm_raster("Value", palette = viridis(n=50), n=50, legend.show=FALSE, legend.hist = TRUE, legend.hist.title = "Species Richness") +
tm_shape(mt) +
tm_borders("white", lwd = .75) +
tm_shape(summary.df.inc) +
tm_polygons(col = "meaninc", border.col = "white", title="Mean Income") +
tm_legend(outside = TRUE)
Task 4: Build a choropleth with ggplot
Your map should have a basemap, should rely on more than one aesthetic (color, transparency, etc), and combine multiple layers.
ggplot(summary.df.pov) +
geom_sf(mapping = aes(fill = family_poverty)) +
geom_sf(data=mt, fill=NA,color="black")
ggplot(summary.df.inc) +
geom_sf(mapping = aes(fill = meaninc)) +
geom_sf(data=mt, fill=NA,color="black")
chloroplath
bg <- ggmap::get_map(as.vector(st_bbox(mt))) # basemap## Source : http://tile.stamen.com/terrain/7/22/43.png## Source : http://tile.stamen.com/terrain/7/23/43.png## Source : http://tile.stamen.com/terrain/7/24/43.png## Source : http://tile.stamen.com/terrain/7/25/43.png## Source : http://tile.stamen.com/terrain/7/26/43.png## Source : http://tile.stamen.com/terrain/7/27/43.png## Source : http://tile.stamen.com/terrain/7/22/44.png## Source : http://tile.stamen.com/terrain/7/23/44.png## Source : http://tile.stamen.com/terrain/7/24/44.png## Source : http://tile.stamen.com/terrain/7/25/44.png## Source : http://tile.stamen.com/terrain/7/26/44.png## Source : http://tile.stamen.com/terrain/7/27/44.png## Source : http://tile.stamen.com/terrain/7/22/45.png## Source : http://tile.stamen.com/terrain/7/23/45.png## Source : http://tile.stamen.com/terrain/7/24/45.png## Source : http://tile.stamen.com/terrain/7/25/45.png## Source : http://tile.stamen.com/terrain/7/26/45.png## Source : http://tile.stamen.com/terrain/7/27/45.png## Source : http://tile.stamen.com/terrain/7/22/46.png## Source : http://tile.stamen.com/terrain/7/23/46.png## Source : http://tile.stamen.com/terrain/7/24/46.png## Source : http://tile.stamen.com/terrain/7/25/46.png## Source : http://tile.stamen.com/terrain/7/26/46.png## Source : http://tile.stamen.com/terrain/7/27/46.pngggmap(bg) +
geom_sf(data = summary.df.pov, mapping = aes(fill = family_poverty), inherit.aes = FALSE) +
geom_sf(data=mt, fill=NA,color="black", inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326))## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
ggmap(bg) +
geom_sf(data = summary.df.pov, mapping = aes(fill = Value,
alpha = (landvalrastmtwy - max(landvalrastmtwy, na.rm=TRUE))/(max(landvalrastmtwy, na.rm=TRUE)-min(landvalrastmtwy, na.rm = TRUE))), inherit.aes = FALSE) +
geom_sf(data=mt, fill=NA,color="black", inherit.aes = FALSE) +
scale_fill_viridis(option="magma")+
coord_sf(crs = st_crs(4326))## Coordinate system already present. Adding new coordinate system, which will replace the existing one.## Warning in max(landvalrastmtwy, na.rm = TRUE): no non-missing arguments to max;
## returning -Inf
## Warning in max(landvalrastmtwy, na.rm = TRUE): no non-missing arguments to max;
## returning -Inf## Warning in min(landvalrastmtwy, na.rm = TRUE): no non-missing arguments to min;
## returning Inf## Warning in max(landvalrastmtwy, na.rm = TRUE): no non-missing arguments to max;
## returning -Inf
## Warning in max(landvalrastmtwy, na.rm = TRUE): no non-missing arguments to max;
## returning -Inf## Warning in min(landvalrastmtwy, na.rm = TRUE): no non-missing arguments to min;
## returning Inf
Task 5: Build 2 cartograms
Follow the examples to build cartograms that display your region of interest based on variables other than area.
mt.pop <- cartogram_cont(mt.income, "pop", itermax = 5)## Mean size error for iteration 1: 4.53233979977975## Mean size error for iteration 2: 3.83147348099061## Mean size error for iteration 3: 3.22771833387299## Mean size error for iteration 4: 2.71033007692304## Mean size error for iteration 5: 2.28592163165744mt.inc <- cartogram_cont(mt.income, "medianincome", itermax = 5)## Mean size error for iteration 1: 1.60285632679## Mean size error for iteration 2: 1.3422147510671## Mean size error for iteration 3: 1.20010655792736## Mean size error for iteration 4: 1.12155703205414## Mean size error for iteration 5: 1.07689808829867mt.pov <- cartogram_cont(mt.poverty, "family_poverty", itermax = 5)## Mean size error for iteration 1: 3.76105120437958## Mean size error for iteration 2: 3.0949361986733## Mean size error for iteration 3: 2.56048558715464## Mean size error for iteration 4: 2.13421537943059## Mean size error for iteration 5: 1.80716554652857# cartogram based on population
tm_shape(mt.pop) + tm_polygons("pop", style = "jenks", title="Population") +
tm_layout(frame = FALSE, legend.position = c("left", "bottom")) 
# cartogram based on median income
tm_shape(mt.inc) + tm_polygons("medianincome", style = "jenks", title="Median Household Income") +
tm_layout(frame = FALSE, legend.position = c("left", "bottom")) 
# cartogram based on family poverty
tm_shape(mt.pov) + tm_polygons("family_poverty", style = "jenks", title = "Family Poverty Rates") +
tm_layout(frame = FALSE, legend.position = c("left", "bottom")) 
Question 3: Reflect on the different maps you’ve made, what do the different visualizations tell you about the data you plotted? Base on the data plotted, Montana has higher species richness in the western part of the state and that income in areas surrounding PAs also appears to be higher for the western part of the state. Poverty rates also seemed to follow that same trend.
Question 4: How might you improve the maps you’ve made? I am not sure if I did the baseplot map incorrectly but I found that one the least interesting. Visualizing the data in that way made it seem like being a protected area was the function of poverty, but it would be more interesting to see that if it is a regional thing or truly something correlated to PA status. I would also change the color to be on a red-orange scale since my mind biases blue towarsd something positive which family poverty is not. With more time I would have also created maps to compare things like HMI instead of just land value.
Question 5: Is a map the best way to evaluate the data you worked with? Why or Why not? I think using maps was pretty appropriate for the purposes of this assignment. I was not as impressed with overlaying the data with the opaque blocks such as the species richness poverty maps. These cartograms were more interesting and straight forward to interpret than some of the earlier maps and overall this assignment is useful for assessing some intiial potential trends but looking at more mathematical relationships and visualizing things that way could provide more meaningful data.